Currently only on github - but will be available on CRAN soon
# remove.packages("brainMapR")
devtools::install_github("baptisteCD/brainMapR")
library(brainMapR)Some brain association maps (for age at MRI) are given with the package. They are stored in the package source folder.
system.file("extdata/07_BWAS_MOA_realPheno", "", package = "brainMapR",mustWork = TRUE)
system.file("extdata/03_BWAS_uncorrected_realPheno", "", package = "brainMapR",mustWork = TRUE)for (pheno in c("Age_MRI","body_mass_index_bmi_f21001_2_0","fluid_intelligence_score_f20016_2_0", "sexuses_datacoding_9_f31_0_0","smoking_status_f20116_2_0", "sexuses_datacoding_9_f31_0_0") ){
for (scenario in c("03_BWAS_uncorrected_realPheno", "09_BWAS_ICVagesex_realPheno", "04_BWAS_5globalPCs_realPheno", "05_BWAS_10globalPCs_realPheno", "06_BWAS_10specificPCs_realPheno", "07_BWAS_MOA_realPheno", "10_BWAS_MOA_multiORM_QC_realPheno", "15_BWAS_MOA_allModa_FE_realPheno")){
identifyClustersBWAS(pathToTheBWASresults =system.file(paste0("extdata/", scenario), "", package = "brainMapR",mustWork = TRUE) , bwasFile = paste0("BWAS_",pheno), outputFolder = "../", signifThreshold = 5e-8)
print(paste(pheno, scenario))
}
}Legend to be used in plots
cols=c(RColorBrewer::brewer.pal(n = 10, name = "RdYlBu")[10:6],RColorBrewer::brewer.pal(n = 10, name = "RdYlBu")[5:1]) # Select palette colours
png("legendbar.png", width = 5, height = 15, units = "cm", res = 400)
par(mar=c(2,4,2,2))
my.colors = colorRampPalette(cols)
z=matrix(1:100,nrow=1)
x=1
y=seq(-0.4,0.4,len=100) # supposing 3 and 2345 are the range of your data
image(x,y,z,col=my.colors(20),axes=FALSE,xlab="",ylab="", main="Correlation")
axis(2)
dev.off()library(knitr)
include_graphics(path = "examplePlots/legendbar.png", dpi = 50)Brain surface plots, made in R using the rgl package. The cortical surface can take a while to render but the plotting is very versatile.
The brain plot functions use the summary files created using identifyClustersBWAS().
Brain plots will silently open the phenotype file (located in UKB_phenotypes15K) to scale association betas into correlation coefficients.
library(viridis)
cols=viridis_pal(option = "C")(7)
for (var in c("Age_MRI","body_mass_index_bmi_f21001_2_0","fluid_intelligence_score_f20016_2_0", "sexuses_datacoding_9_f31_0_0","smoking_status_f20116_2_0", "sexuses_datacoding_9_f31_0_0") ){
for (scenario in c("03_BWAS_uncorrected_realPheno", "09_BWAS_ICVagesex_realPheno", "04_BWAS_5globalPCs_realPheno", "05_BWAS_10globalPCs_realPheno", "06_BWAS_10specificPCs_realPheno", "07_BWAS_MOA_realPheno", "10_BWAS_MOA_multiORM_QC_realPheno", "15_BWAS_MOA_allModa_FE_realPheno")){
plotCortical(inputPath = "../", bwasFile = paste0("BWAS_",pheno), pathPhenotypeFile = "../../../UKB_phenotypes15K/Age_MRI.phen", outputPath = " ../")
plotSubcortical(inputPath = "../", bwasFile = paste0("BWAS_",pheno), pathPhenotypeFile = "../../../UKB_phenotypes15K/Age_MRI.phen", outputPath = " ../")
plotSubcortical_flat(inputPath = "../", bwasFile = paste0("BWAS_",pheno), pathPhenotypeFile = "../../../UKB_phenotypes15K/Age_MRI.phen", outputPath = " ../")
combineCorticalSubcorticalPlots(inputPath = "../", bwasFile = paste0("BWAS_",pheno), pathPhenotypeFile = "../../../UKB_phenotypes15K/Age_MRI.phen", outputPath = " ../" )
}
}Left cortical surface area associations with age - GLM without covariates
library(knitr)
hemi="lh"
moda="area"
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_inside.png"))
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_outside.png")) > Right cortical surface area associations with age - GLM without covariates
library(knitr)
hemi="rh"
moda="area"
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_inside.png"))
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_outside.png"))Left cortical thickness associations with age - GLM without covariates
library(knitr)
hemi="lh"
moda="thickness"
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_inside.png"))
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_outside.png"))Right cortical thickness associations with age - GLM without covariates
library(knitr)
hemi="rh"
moda="thickness"
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_inside.png"))
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_outside.png"))Sucbcortical surface area associations with age - GLM without covariates
library(knitr)
hemi="rh"
moda="LogJacs"
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_inside.png"))
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_outside.png"))
hemi="lh"
moda="LogJacs"
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_inside.png"))
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_outside.png"))include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/Plots_CombinedAge_MRI.png"))library(viridis)
cols=viridis_pal(option = "C")(7)
BrainMapManhattanPlot(inputPath = system.file("extdata/07_BWAS_MOA_realPheno", "", package = "brainMapR",mustWork = TRUE), path = "03_BWAS_uncorrected_realPheno/", bwasFile = "BWAS_Age_MRI.linear", yMax = 120, signifThreshold =5e-8, phenotypeLabel = "Age")include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/Manhathan_BWAS_Age_MRI.linear_simple.png"))UKB_BWASvars_labels.txt contains 3 columns: variable names, category (e.g. cognition, lifestyle) and variable label (used for plots and tables)
pheno=read.table("UKB_BWASvars_labels.txt", header=T, stringsAsFactors = F)
pheno=pheno[which(pheno$variable %in% c( "body_mass_index_bmi_f21001_2_0","fluid_intelligence_score_f20016_2_0", "Age_MRI","sexuses_datacoding_9_f31_0_0", "smoking_status_f20116_2_0")),]
vall=NULL
for (iii in 1:5){
all=NULL
for (scenario in c("03_BWAS_uncorrected_realPheno","09_BWAS_ICVagesex_realPheno", "05_BWAS_10globalPCs_realPheno", "07_BWAS_MOA_realPheno", "15_BWAS_MOA_allModa_FE_realPheno", "10_BWAS_MOA_multiORM_QC_realPheno")){
vres=read.csv(paste0("../../",scenario, "/", "BWAS_signif_", pheno$variable_clean[iii], ".csv" ))
res=read.table(paste0("../../",scenario, "/", "Results_clusterFWER_", pheno$variable[iii], ".txt" ), header=T)
nclust=res[,c("NumberClusters_thickness", "NumberClusters_area", "NumberClusters_thick", "NumberClusters_LogJacs")]
maxclust=res[,c("maxSizeCluster_thickness", "maxSizeCluster_area", "maxSizeCluster_thick", "maxSizeCluster_LogJacs")]
all=rbind(all, dim(vres)[1], sum(nclust), max(maxclust, na.rm = T))
}
vall=cbind(vall, all)
}
colnames(vall)=pheno$variable_clean
writexl::write_xlsx(as.data.frame(vall), "../../Summary_BWAS_realPheno_bonferroniNtraits_R1IEEE.xlsx")
# N clusters cortex
vall=NULL
for (iii in c(1,2,3,5)){
all=NULL
for (scenario in c("03_BWAS_uncorrected_realPheno","09_BWAS_ICVagesex_realPheno", "05_BWAS_10globalPCs_realPheno", "07_BWAS_MOA_realPheno","15_BWAS_MOA_allModa_FE_realPheno", "10_BWAS_MOA_multiORM_QC_realPheno")){
#vres=read.csv(paste0("../../",scenario, "/", "BWAS_signif_", pheno$variable_clean[iii], ".csv" ))
res=read.table(paste0("../../",scenario, "/", "Results_clusterFWER_", pheno$variable[iii], ".txt" ), header=T)
nclust=res[,c("NumberClusters_thickness", "NumberClusters_area", "NumberClusters_thick", "NumberClusters_LogJacs")]
maxclust=res[,c("maxSizeCluster_thickness", "maxSizeCluster_area", "maxSizeCluster_thick", "maxSizeCluster_LogJacs")]
all=rbind(all, sum(nclust[,1:2]))
}
vall=cbind(vall, all)
}
colnames(vall)=pheno$variable_clean# Brain plot
source("RR_9.1_BWAS_manhathanPlot_otherFunctions.R")
library(Morpho)
library(plyr)
library(viridis)
for (scenario in c( "07_BWAS_MOA_realPheno")){
# Plot cortical and subcortical
plotCortical_noSignif(path = paste0("path/to/working/directory/", scenario, "/"), variable = "body_mass_index_bmi_f21001_2_0")
plotSubcortical_noSignif(path =paste0("path/to/working/directory/", scenario, "/"), variable = "body_mass_index_bmi_f21001_2_0")
combineCorticalSubcorticalPlots_noSignif(path = paste0("path/to/working/directory/", scenario, "/"), variable = "body_mass_index_bmi_f21001_2_0")
}A work by by [Baptiste Couvy-Duchesne] - 07 October 2021